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ABSTRACT 

Generally, wave field reconstructions obtained by phase- 
retrieval algorithms are noisy, blurred and corrupted by var- 
ious artifacts such as irregular waves, spots, etc. These dis- 
turbances, arising due to many factors such as non-idealities 
of optical system (misalignment, focusing errors), dust on 
optical elements, reflections, vibration, are hard to be local- 
ized and specified. It is assumed that there is a generalized 
pupil function at the object plane which describes aberra- 
tions in the coherent imaging system manifested at the sensor 
plane. Here we propose a novel two steps phase-retrieval al- 
gorithm to compensate these distortions. We first estimate the 
cumulative disturbance, called "background", using special 
calibration experiments. Then, we use this background for 
reconstruction of the object amplitude and phase. The second 
part of the algorithm is based on the maximum likelihood 
approach and, in this way, targeted on the optimal amplitude 
and phase reconstruction from noisy data. Numerical ex- 
periments demonstrate that the developed algorithm enables 
the compensation of various typical distortions of the optical 
track so sharp object imaging for a binary test-chart can be 
achieved. 

Index Terms — Noise in imaging systems, Spatial light 
modulators, Phase retrieval, Inverse problems 

1. INTRODUCTION 

The phase contains important information on the shape of the 
object, which is useful in metrology and 3D imaging, e.g. mi- 
croscopy, astronomy, material analysis, etc. The conventional 
sensors detect only the intensity of the light. Since the phase 
cannot be measured directly and it is systematically lost in 
observations, computational phase recovering techniques are 
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required for imaging and data processing. Phase recovering 
and, in general, the reconstruction of the object amplitude and 
phase is referred to as the phase-retrieval problem. 

Perhaps the first iterative method for phase retrieval from 
intensity measurements were the well-known Gerchberg- 
Saxton algorithm [1], initially for a single observation plane, 
and its variation devised by Misell [2] for two defocusing 
images at different measurement planes. The idea consisting 
in the iterative replacement of the estimated magnitude by 
measured and prior information is further developed for var- 
ious applications by many authors (e.g. [3, 4, 5, 6]). Similar 
methods are proposed for Fresnel instead of Fourier trans- 
forms as the transfer functions of the wave field propagation 
[7, 8, 9]. Various phase-retrieval algorithms based on these 
landmark works are systematized by Fienup [10] introducing 
classical types of iterative phase-retrieval algorithms. Mul- 
tiple measurements gain an observation redundancy that can 
be exploited in order to improve the quality of the complex- 
valued object reconstruction [11, 12, 13]. 

The above imaging techniques are mainly based on an 
ideal wave field propagation modeling derived from the scalar 
diffraction theory [14]. In practice, wave fields in real coher- 
ent imaging systems and their observations are quite different 
from those predicted by theory, hence wave field reconstruc- 
tions obtained by simulations (i.e. theoretical results) and 
using real experimental data can dramatically vary. The re- 
constructions obtained from the real data differ from simu- 
lated ones by multiple and well seen artifacts which can have 
a form of disturbed background with irregular waves, spots, 
random noise, etc. These systematic distortions appear due 
to many factors such as non-idealities of optical system (mis- 
alignment, focusing errors, aberrations), dust on optical ele- 
ments, reflections, vibration, etc. 

In this paper we consider a 4f optical system with an SLM 
located across the Fourier domain of the first lens. This sys- 
tem is used for capturing multiple intensity observations at the 
sensor plane for the phase reconstruction (see Fig. 1, [15]). 
The reconstruction from this data is very sensitive to all dis- 
turbing factors because it is an ill-posed problem. One of the 




sensor 



Fig. 1. Experimental setup of the 4f optical system used for record- 
ing measurement data [15]. The lenses L\ and L2 in the 4f con- 
figuration provides an accurate mapping of the object wave field to 
the parallel observation (sensor) plane. An optical mask with the 
complex-valued transmittance M r located at the Fourier plane (a 
phase modulating SLM) enables linear filter operations. 



strongest sources of disturbances is the used SLM due to a 
high sensitivity of 4f system to modulation of the wave field 
at the Fourier plane. 

In general, there are diversity of numerical approaches, 
which are used for calibration [16], filtering parasitic reflec- 
tions [17], compensating for curvature introduced by micro- 
scope objective [18, 19], for aberrations [20] or astigmatism 
[21]. In this work we follow an essentially different idea. It is 
assumed that there is a generalized pupil function [14, §6.4.1] 
at the object plane which describes aberrations in the coher- 
ent imaging system manifested at the sensor plane. Namely, 
we do not try to identify particular sources of the disturbances 
but estimate and compensate their accumulated effects by re- 
calculating them to the entrance pupil of the used 4f config- 
uration. In the following, the cumulative distortions are re- 
ferred to as "background" distortions. Thus, we are devel- 
oping two step phase-retrieval algorithm: firstly, we estimate 
this background disturbance using special calibration experi- 
ments and then use it to reconstruct the object amplitude and 
phase. In this work we apply the variational constrained max- 
imum likelihood formulation with parallel processing of mul- 
tiple intensity observations proposed in our previous works 
[22, 23, 24, 25]. Moreover, we incorporate prior information 
on the true object wave field: in our experiments we recon- 
struct a binary object with unknown lower and upper levels. 

Let uq(x), x <E R 2 be a true object wave field at the 
entrance pupil of the system. Taking into consideration the 
non-ideality of the optical system, we introduce a "disturbed" 
object wave field u (x) as a product of a typically unknown 
background (cumulative distortion) wave field ub by the true 



object wave field uq(x) as 

u (x) =u (x) -u B (x), (1) 

where the diacritic " emphasizes the difference of the cor- 
rupted wave field uo from the ideal one uq. The standard 
phase-retrieval techniques are able to give the reconstruction 
of the disturbed wave field uq (x) only and they are not able 
to separate the background in order to estimate the true wave 
field u (x). 

In our work, we try to reconstruct the disturbances by per- 
forming a calibration procedure (to estimate the background 
ub) and use it to extract the true object uq. At first glance, 
this problem looks trivial: one may produce the experiments 
with a known invariant uq(x), for instance uq(x) = 1, ob- 
tain the estimate ub (x) and then recalculate the estimate for 
the object as iio(x) = uq(x)/ub(x). However, a priori infor- 
mation about the object in the used sparse modeling concerns 
the true object wave field uo but not the disturbed one uq. 
Thus, we are processing the recalculated object estimate at 
each iteration (uq(x) = u t Q (x)/uB(x), t = 0,1,2,...), and 
the structure of the developed iterative phase-retrieval algo- 
rithm is therefore essentially different from the trivial guess. 

The paper is organized as follows. In Section 2, the image 
formation in a 4f optical system and the observation model 
are presented. The constrained variational approach for the 
phase retrieval and the sparse modeling for the object phase 
and amplitude are introduced in Section 3. The proposed 
phase-retrieval algorithm with the background compensation 
is presented in Section 4. Numerical experiments for the ob- 
ject wave field reconstructions from real data are shown and 
discussed in details in Section 5. 

2. OBSERVATION MODEL 

Let us consider the image formation model in a conventional 
4f configuration of the coherent imaging system linking com- 
plex amplitudes at the object and measurement planes. 

Let us denote complex amplitudes at the object and mea- 
surement (sensor) planes by uq(x) and u r (y), respectively. 
The lenses L\ and L2 with the focal length / arranged in the 
4f configuration provides an accurate mapping of the object 
wave field into the parallel measurement plane. A reflective 
phase modulating spatial light modulator (SLM) is placed at 
the Fourier plane of the first lens [15]. The used 4f optical 
system is illustrated in Fig. 1 . 

Let us assume for a moment that there are no distortions 
in the optical track. It is well known that the link between the 
wave fields at the object uq (x) and the Fourier planes uf ( ) 
is given as follows [14] 

where J-{-} denotes the 2D integral Fourier transform, A is a 
wavelength. 




If the optical mask (SLM) inserted at the Fourier plane has 
the complex- valued transmittance M r (jj), then the output of 
the optical system is defined as 



Similarly, discrete model for Eq.(3) has the form 



u r (y) 
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(3) 



All these wave field distributions are given in the 2D lat- 
eral coordinates: here we use the variables x, y, v G M 2 for 
the object, sensor and Fourier planes, respectively. 

2.1. Discrete modeling 

For discrete modeling, the continuous arguments are replaced 
by the digital ones with a corresponding replacement of the 
continuous functions by their discrete counterparts: Uo(x) — > 
u a (k 1 A Xl ,k 2 A X2 ), u r (y) -» u r (hA yi , h\ 2 ), u F {fj) -> 

Bfl^rjj, ^fV 2 ) w i m 2-D integer arguments k — fc 2 ), 
I = (h,h) an d V — (^i?^)- This discretization is dictated 
by the use of a digital camera and a pixelated SLM as a 2D ar- 
ray of liquid crystal cells. Thus, we hereafter consider the dis- 
crete wave fields at the object uo[ki, k 2 ], Fourier up [r] 1 ,ri 2 ] 
and sensor planes u r [li,l 2 ] with various pixel sizes x 



where F{-} in Eqs. (5) and (6) denotes the 2D DFT op- 
erator, o stands for the Hadamard (elementwise) product 
and M r [rj 1 , rj 2 ] is the discretized optical mask at the Fourier 
plane. Note that the wave field propagation can be easily 
realized much faster using FFT (cf. [24, 25]). 

Taking into account the vector-matrix notation and the 
distortions in the real optical system (Eq. (1)), the forward 
wave field propagation from the object to the sensor plane 
can be given in the form 



A r ■ u , r = 1, ...K, 



(7) 



where u € C" xl is a complex- valued vector, correspond- 
ing to the disturbed object discrete 2D wave field distribution, 



= N x ■ N v . A r g 



is a forward propagation opera- 



tor corresponding to the optical mask M r at the Fourier plane 
(programmed using the SLM), and if is a number of these 
various optical masks. 



A Xl , A Vl x A Vl and A yi x A yi , respectively. In general, 

these images can be rectangular of different size N Xl x N Xl , 2 .2. Noisy intensity observations 



N Vl x N V2 and N yi x N yi , respectively 

We use a vector-matrix notation for complex-valued dis- 
tributions of the wave fields. 2D discrete distributions (ma- 
trices) are vectorized to the complex-valued column vector 
[26]. Bold lower case characters are used for the vectors. 
Matrices are defined by bold upper case to distinguish them 
from vectors. Thus, Uo[fc], u F [r]] and u r [l] are column vec- 
tors constructed by vectorization of the corresponding 2D dis- 
crete wave field distributions at the object Uo[fci, k 2 ], Fourier 
Uf^, 7? 2 ] and sensor planes XJ r [h, h], respectively. 

In this work we assume that the pixel size at the object 
and sensor planes is the same (A Xl — A yi , A X2 — A y2 ) and 
images Uo, Uj? and U r are of the same size N x x N y for all 
planes. Let us also assume that the following conditions are 
fulfilled [27] 



Assume that we have a set of K experiments produced with 
different masks {M r },., r = 1,...,K. The problem is to 
reconstruct a complex-valued true object wave field u from 
multiple noisy intensity observations {o r } measured at the 
sensor plane. These measurements are represented in the 
vector-matrix notation following to Eq. (7) as follows 



o r [Z] = |u r [Z]| 2 + £,-. 



1,...K, 



(8) 



where the noise is assumed to be zero-mean Gaussian with 
the variance of, e r [l] ~ Af(0, cr 2 ), independent for different 
I and r. The observation vectors {o r } correspond to the 2D 
distributions on the regular discrete grid located at the sensor 
plane. 



A V1 A X1 N X 
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Then, discretization of the integral in Eq. (2) defines 
U_f [Vi j V2] as 2D discrete Fourier transform (DFT) of Uo [fci , k 2 ] 
in the form 



(5) 



A X1 A X2 



fcl = ft 



E 



\J (kiA xl ,k 2 A X2 ) x 



-k 2 = - 



x exp( T 



2- 



iXf 



-N„ 



< 



(rj 1 A Vl k 1 A Xl 



■q 2 A V2 k 2 A X2 )) 



< 



*vi 

2 



1. 



N„ 



< 



'72 



< 



3. SPARSE OBJECT MODELING AND 
VARIATIONAL FORMULATION 

It is assumed in sparse modeling approach that the "true" ob- 
ject distribution Uo can be approximated by a small number 
of non-zero elements of basis functions. The ideal basis func- 
tions for the object approximation are unknown a priori and 
selected from a given set of potential bases (dictionaries). In 
general, we deal with a complex-valued object wave field, 
and hence consider nonlinear modeling of the wave field with 
separate approximations for the object phase and amplitude 
[22, 23, 24, 25]. We represent the object wave field in the 
form u = a o exp(j • tp ), where a = abs(u ) e M" 
and (p = angle(uo) £ K" denote the object amplitude and 



phase, respectively. Sparse object approximation can be given 
in the analysis or synthesis form as follows 



6 a = * a • abs(u ),9 v = 
a = * Q ■ O a ,ip 



angle(\io) 

■ din 



(analysis) 
(synthesis) 



(9) 

Here ^ a ,^ v and 3> a , 3?^ are the frame transform matri- 
ces, and the vector a ,9 v £ R m can be considered as a spec- 
trum (m ^> n) in a parametric data adaptive approximation. 
Subindices a and 93 are shown for the amplitude and phase, 
respectively. It is recognized that, in contrast to classical or- 
thonormal bases (to = n), overcomplete frame based model- 
ing is a much more efficient for imaging [28, 29] and results 
in a better wave field reconstruction accuracy. The sparsity 
of approximation is characterized by either the £q norm |o 
defined as a number of non-zero components of the vector 
or the £1 norm as a sum of absolute values of components 
of the vector = Yl s \® s \- ^ sma U er value of the norm 

means a higher sparsity of approximation. Note that results 
obtained by £0 or l\ norms are shown to be closed to each 
other [30], what allows replacing the nonconvex £ norm by 
the convex l\ norm in many variational settings. 

The main intention is to find sparsest (shortest) models 
for phase and amplitude with smallest values of the £q or £\ 
norms. The separate sparse modeling for the object phase 
and amplitude is realized via the powerful BM3D-frame filter, 
specified for denoising and other imaging problems [31, 32, 
33]. 

Assume that the background wave field ub is given. Tak- 
ing into account the sparse modeling for the object amplitude 
and phase, the wave field reconstruction is performed by min- 
imization of the criterion J 



r— 1 



■ 2 1 1 2 



•Ta||0a||p+<7v||0J|p (10) 



subject to u r = A r • iio, r = 1, ...K, (11) 

u = u o u B (12) 

a = * a • a&s(u ), 8 V = & v • angle(uo), (13) 

a o = V a d a ,tp = i& v ,e <p (14) 

where 1 1 ■ 1 1 2 stands for the Euclidean norm and regularization 
terms for phase and amplitude are taken using the l p norms 
(p = {0, 1}). The positive parameters r a and t v in Eq. 
(10) define a balance between the fit of observations, smooth- 
ness of the wave field reconstruction and the complexity of the 
used model (cardinality of spectra 9 a , 6 V of the object am- 
plitude and phase). Note that the constraint for the forward 
wave field propagation (1 1) is presented for the disturbed ob- 
ject wave field iio (Eq. (12)), and the used sparse modeling is 
given for the true object: namely, the analysis (13) and syn- 
thesis (14) are calculated using the frame transform matrices 
for the compensated object u [k] = Uo[fc]/u#[fc] (rather for 
the object amplitude and phase). 



3.1. Multi-objective optimization 

It is shown in [33] that a multi-objective optimization can be 
much more efficient than the minimization of the single crite- 
rion J due to a simpler implementation (filtering and inverse 
procedure are decoupled) and resulting better reconstruction 
quality. Thus, instead of the constrained minimization of (10) 
we arrive at the unconstrained minimization of two criterion 
functions J\ and J 2 with changing the constraints for sparse 
modeling by the quadratic penalties with positive weights 



Jx = — ||uo-v ||i+ (15) 



7o 



K 



•Si^llOr-lV"" ' 

r=l 



2 \\l + — |K - A r -Uol^ 



7r 



Re{Af • (u r - A. r • u )}, 



Jl = T a \\0a\\p + ^ ||6> a - * Q • abs(u )\\ 2 2 + (16) 

+tJ\0<p\\ p + x—\\0 v - * v • ang/e(u )|||, 



where (-) H stands in Eq. (15) for the Hermitian conju- 
gate, Vo = ^a^a ° exp(j ■ tyipO^p) is an approximation of 
the complex- valued object distribution Uo, vo = vo o ug. 
{A r } £ C n are the complex-valued vectors of the Lagrange 
multipliers (see [34]). Note that the linear and quadratic 
penalties related to the forward propagation are involved with 
the same positive parameters j r . The analysis and synthesis 
constraints in Eqs. (13) and (14) are replaced by the penalties 
with the corresponding positive parameters r ) a , j v and 7 , 
in Eqs. (15) and (16), what is a standard tools to deal with 
constrained optimization [35]. 

Note that the criterion function J 2 is separable with re- 
spect to a and 8^, thus it can be rewritten as J72 = J2,a + 
J 2 ^, where 



J2,a{ 6 a,abs{\i )) 
1 , 



(17) 



= r -Ileal 



27 tt 



6 a - $ a • abs{u )\\l 



J 2 , v {9 v ,angle(u )) 



(18) 



• \\0 



<p\\p 



angle(u )\\l 



In contrast to [24], where only a quadratic penalization is 
used for the forward propagation model u r = A r -iio, here we 
use the linear Lagrangian summand Rc{A^ • (u r — A r • iio)} 
in Eq.(15). Thus, J\ becomes the augmented Lagrangian ob- 
jective function at least with respect to this particular con- 
straint. 



4. PROPOSED ALGORITHM 

Following our assumption that the object wave field is de- 
graded by the distortions accumulated in the background 
wave field ub (uo = Uo o cf. Eq. (1)), special calibration 
experiments are produced in order to estimate these distur- 
bances. Thus, it is assumed that two sets of experiments are 
made consistently under the same conditions. These experi- 
ments result in two sets of observations {of} and {o r } used 
to estimate the background u# and the the object u , respec- 
tively. The flowchart of the proposed two steps algorithm is 
shown in Fig. 2. 



of. The estimate ub can be computed using the following 
algorithm 

Repeat fort = 0,1, 2,... (20) 

Repeat for r = 1, ...K 

u* = arg min J AL (of , u' B , u r , A' , ) 

Al +1 = A* + a r ■ (u* - A r • u' B ) 
End on r 

u^ 1 = argmin l 7 Ai ({of }, u B , {u*}, {A*}), 

us 

End on t 



4.1. Background reconstruction 

Firstly, we record a number of intensity observations {of} 
corresponding to various optical masks {M, } for the free 
space object (test-image Uo[A;] = 1) and find a complex- 
valued estimate ub (uq = u B ) by optimization of the cri- 
terion function 

K 1 1 

J AL = (x\\n B \\ 2 2 + £ -3 Mo? - |u,| 2 ||^ + (19) 

r—1 

1 2 

H ||ur - A r UE>||2 H Rc{Af (u r - A r ue)} 

Tr it 

as it is described in [34] . The main difference of the structure 
of Jal from J\ consists in the last quadratic penalty term, 
/i in Eq. (19) is the Tikhonov regularization parameter which 
defines a balance between the prior information on ub and the 
fitting of calculated intensities |u r | 2 to the given observations 




u =u °u B o r =\A r u Q \ 2 +£ r 



Fig. 2. A flowchart of the two steps phase retrieval technique 
with the background compensation. The upper block high- 
lighted by a dashed line represents the background calibration 
procedure, where the complex-valued estimate of u B is found 
by AL [34]. The reconstruction of the object using the back- 
ground estimate is obtained by the proposed SPAR — BC 
method. 



This algorithm without background compensation and ad- 
ditional object filtering is identical to the AL algorithm orig- 
inated in [34] but presented with respect to the background 

Vi B . 

This stage is shown in the upper block of Fig. 2 high- 
lighted by a dashed line. 

4.2. Object reconstruction 

Secondly, we record intensity measurements for an object 
{o r } using the same optical masks {M r } as before and re- 
construct the true object wave field using the found back- 
ground estimate u b ■ 

According to the general idea of the multi-objective opti- 
mization, where the alternating minimization of J\ with re- 
spect to Uo = Uo o ub, {u r } and minimization of J2 with 
respect to 6 a , 6 V are used [32, 33, 24], we arrive at the fol- 
lowing iterative algorithm 

u* [fe] = u [fc]/u B [fc] (21) 
61 = argmin J2 :a (0 a ,a6s(u o )) (22) 

a 

6^ = a,rgmmj 2 , v (0 v ,angle(u t o )) (23) 

v* = * o 0* o exp(j • * v 0%) o u B (24) 
u' = argmin Ji(o r ,u ,u r , A',Vq), (25) 

11, 

A* +1 = A* + a r • (u* - A r • u ), r = 1, ...K, (26) 
u +1 = argmmJi({o r },Uo,{u*},{A*},v ) (27) 

Here we first estimate the true object (Eq.(21)) by com- 
pensation of the disturbed u with the background estimate 
ub- It results in the object amplitude and phase estimates. 
Then, Eqs. (22) and (23) enable the spectrum estimates of the 
object amplitude and phase by thresholding in the BM3D- 
frame domain with the thresholds r a r y a and r v 7 v , respec- 
tively [31, 32]. Eq. (24) corresponds to the synthesis of the 
approximation v of the disturbed object from the calculated 
spectra for the true object amplitude and phase and using the 
background u B . Together the operations in Eqs. (22)-(24) 
related to the optimization of J2 can be rewritten in a more 



compact form as follows 

&l +1/2 = BAttD a (abs(ul)), (28) 
<p f +1/2 = BMZD v {angle(vfc)), 

-t t+l/2 , ■ t+l/2x 

v = a o exp(j ■ <p ) o u B 

where BMiD{ ) denotes hereafter the processing by the 
BM3D filter, and the corresponding subindices a and cp em- 
phasize that the filtering is performed with different param- 
eters and different transform matrices \&. and 4>. for the 
amplitude and phase, respectively. In our implementation the 
analysis and synthesis operations, the thresholding and cal- 
culation of the matrices \&. and <fr. are integrated in a single 
block which is called BM3D filter. 

Eqs. (25), (27) are the optimization steps for J\\ the com- 
putation of the complex- valued wave field estimates {u r } at 
the sensor planes and the disturbed object iio from noisy in- 
tensity observations {o r }. The update of the Lagrange vari- 
ables A r is shown in Eq. (26). This second stage is illustrated 
in Fig. 2 under the mentioned block for the background esti- 
mation. 

4.3. Sparse Phase Amplitude Reconstruction with Back- 
ground Compensation {SPAR — BC) 

Following the mentioned two main steps - Eqs. (20) for 
the background estimation and Eqs. (21)— (27) for the 
true complex-valued object extraction - we formulate the ad- 
vanced phase-retrieval approach with background compensa- 
tion. Taking into account Eqs. (28), the reconstruction of 
the true object wave field is performed by the proposed itera- 
tive algorithm called Sparse Phase Amplitude Reconstruction 
with Background Compensation (SPAR — BC). 

The initialization concerns the calculation of the back- 
ground (according to Eqs. (20)), the initial guess for the 
disturbed object Uq = Ug m * (say, again by AL [34]) and 
Lagrange multipliers (e.g. A°[fc] = 0). Note that the trans- 
form matrices for both the synthesis ^ a ^ v and analysis 
&a,&.p may be constructed only once during the initializa- 
tion procedure or may be periodically updated. In this work 
we calculate these matrices only ones for the compensated 
object amplitude abs(\io[k]) / abs(uB[k}) and phase estimates 
angle(u.Q[k]) — angle(u.B[k]). Note also that in our experi- 
ments we use the ^i-norm in the sparse object approximation 
("soft" thresholding of the used BM3D filter, see [25, 33]). 

Note that the output of the SPAR — BC phase-retrieval 
algorithm is not the estimate of the disturbed iio (Step 7), but 
the estimate of the true object wave field Uo calculated in Step 
1 . The derivations of main steps of SPAR — BC (the min- 
imization in Eqs. (22), (23), (25) and (27)) can be found in 
our previous works [34, 24, 25]. Step 4 returns the wave field 
u* +1 ^ 2 at the sensor plane corresponding to the forward prop- 
agation model with the optical mask M r . Step 5 gives the 
updates of u* +1 ^ 2 by their fitting to the observations o r . The 



Algorithm: SPAR - BC 
Input: {of }f =1 , {o r }f =1 
Initialization: u B , Uq, {A"} 
Repeat for t = 0, 1,2... 
1. Object update (background compensation): 
u£[fc]=u*[fc]/u B [fc] 
2. BM3D filtering: 
a* +1/2 = BM3D (a6s(i4)), 

<£>o +1/2 = BM3D v (angle(u$)) 
3. Object approximation synthesis: 

t+l t+l/2 / . t+l/2, 

V = a o °exp(j-y ) 
Repeat for r = 1, ...K 

4. Forward propagation: 
u I t - +1/2 = A,. ■ ul 

5. Fitting to observations: 

u^[l] =g(o r [l],ut +1/2 [l], A* [I]) VI 

6. Lagrange multipliers update: 
A* +1 =A t r + a r - (u* +1 - u* +1/2 ) 
End on r 
7. Disturbed object update: 
f,*+i- fV K -J— A ff A +-L-T )- 1 x 

u — \2-ir=l 7 a 2 A r A r T ^ o i-nxn) x 

x Ef=i • K+ 1 + A* ) + i • (u B o v^ 1 ) 

End on t 



operator defining this update is denoted as Q and described in 
[34, Appendix A]. 

It is shown in [24] that the object reconstruction with 
BM3D filtering can be realized without Lagrange multipliers. 
However, it is found (see [25]) that {A r } help to recover 
small details of the object. In Step 6 {A* } are updated with 
the step a r , and in our experiments we take a fixed step 
a r — a = 1/20 for all K observations. 

In this work we use the same noise variation at all sen- 
sor planes (a 2 , = a 2 ) and take the equal parameters for the 
Lagrangian multipliers, j r = 7. Then, it is easy to see that 
the estimate Uq computed in Step 7 of SPAR — BC consists 
of two parts: the disturbed object estimate calculated from the 
observations and the filtered object approximation found from 
the output of the BM3D filter. Then, this step of the algorithm 
can be given in the form 

K 

u t + 1 = ^B r .(u<+ 1 +A<) + / t -vS+\ (29) 

r=l 

where Vq +1 = u_b v o +1 an d me transform matrix B r is 
given in the form 

K 

B. r = (^AfA r + «-I„ xn )- 1 Af, (30) 

r=l 

and k = 7<7 2 /j . In particular, for all our experiments k is 
equal to 90/25. 



5. NUMERICAL RESULTS 

In this Section a high performance of the proposed algorithm 
is demonstrated by example of amplitude reconstructions 
from real experimental data. 

5.1. Binary object model 

Here we consider reconstruction of a binary object with the 
amplitude given as 



ao[fc] = abs(uo[k] 



p x ,foxk elid, 
/3 , for k £ X\X X , 



(31) 



where X is a support of the image, f3 £ R + and f3 1 £ R+ 
stand for the lower and upper level of the object amplitude 
signal, respectively. The set X\ defines the indices of the up- 
per level and the set X = X\Xi defines the indices of the 
lower level. Both the levels j3 , j3 1 and the sets X , Xi are 
unknown and should be reconstructed. The U.S. Air Force 
resolution test-chart is used for ao. For the amplitude-only 
object the phase should be equal to zero, angle(uo[k]) = 0. 
In practice, the laser beam passing through the chart under- 
goes some phase transformations. These transformations de- 
fine the phase characteristics of the object Uo, which are un- 
known. The only thing which can be stated is that the pattern 
of the object phase reflects the binary features of the ampli- 
tude model (31). Thus, we are looking for the complex-valued 
object u . 

5.2. Settings of parameters 

The standard settings of the phase-retrieval problem assume 
a multi-plane lensless system with varying distances between 
the parallel object and sensor planes. The intensity measure- 
ments at the sensor planes are used for reconstruction of a 3D 
wave field including both the phase and amplitude. Follow- 
ing [15] the considered 4f optical system works as an imitator 
of this multi-plane lensless scenario. The principal difference 
is that the sensor plane is immobile and fixed at the distance 
4/ from the object plane. The effect of the varying distances 
is obtained by a phase modulating SLM located at the Fourier 
plane, where different optical masks M r corresponding to the 
propagation distances are programmed as 



(32) 



. 7T 

exp(2z— z r 
A 




In Eq. (32) z r = Z\ + (r - 1) • A z , r = 1, ...K are 
the distances between the object and sensor planes. In our 
experiments K=5, Zi=20mm is the distance from the object 
to the first measurement plane, A z =2mm is the fixed distance 
between successive measurement planes. 



Due to the bandlimitedness (fixed size of the SLM) and 
discrete representation of the optical mask (on a 2D array of 
liquid crystal cells of the SLM), the experimental results at 
the sensor plane will be different from the output of the model 
calculated using the angular spectrum decomposition (ASD, 
[14]). These differences are consider as components of the 
background to be estimated and compensated. 

In our discrete wave field propagation model, the pixels at 
the sensor and Fourier planes are square of the different size 
A yi x A yi =3.45x3.45 (pm) and A Vl x A„ 2 =8x8 (/xm), re- 
spectively, with 100% fill factor [36, 37]. The object is pixe- 
lated with the sensor size pixels: A Xl = A X2 = A yi = A V2 . 
The transparent U.S. Air Force resolution test-chart (MIL- 
STD-150A) inserted in the front focal plane of the first lens 
L\ is illuminated by collimated coherent light with wave- 
length A=532 nm (i.e. a green Nd:YAG laser is used). The 
employed SLM was supplied by Holoeye Photonics AG and 
configured to provide full 2ir phase modulation. The focal 
distance of lenses used in the 4f configuration is /=150 mm 
what equates with the image size 2892 x 2892 pixels accord- 
ing to the sampling conditions, Eq. (4). The measurement 
area is smaller and here we reconstruct only a part of the ob- 
ject of the size N x x N y (2048 x 2048) for the corresponding 
computational focal distance / c =106.25 mm (see Eqs. (5)- 
(6)). Note that /=150 mm in Eq. (32) defining the optical 
masks M r . 

The algorithm is implemented for a graphic processing 
unit (GPU) in order to use the advantage of parallel process- 
ing of u' +1 and Uq +1 . The GPU realization results in a sig- 
nificant acceleration what is crucial especially for large im- 
ages [25]. The presented results are computed in MATLAB 
7.13 (R2011b) using GPU Nvidia GF460GTX with CUDA 
4.1. The computer used for experiments is Intel i5 2500 (4 
physical cores) at 3.3 GHz; 8Gb RAM, Windows 7 SP1. 

5.3. Modification of SPAR BC for binary object 

A special modification of the filtering procedure is developed 
targeted to improve reconstruction of a binary object. The 
BM3D filtering (Step 2 of SPAR - BC) is replaced by 

a^, +1/3 = BM3D a (abs(ul) - ft) + ft (33) 
a* +1/2 = BM3A,(«4 +1/3 -#)+# 

where (3 and j3\ are scalar variables. These variables are 
calculated as medians of ag = a&s(u ) over the sets: Xq = 
H : < a < p*} andX{ = {a* : a > p*} 



median a t eX t(ao), 



(34) 



In order to estimate these classes Xq and X\, correspond- 
ing to small and large values of a^, we use the thresholding 
parameter p l calculated using the Otsu algorithm [38]. In the 



procedures (33), successive subtractions of f3 and /? x makes 
the image flatter first in the area of low values of binary ampli- 
tude signal and after that in the area of its high values. Exper- 
iments show that this flattening enables much more efficient 
filtering of artifacts for the estimate of ao in case of binary 
object. 

For the phase filtering we make the flattering procedure 
simpler because for the considered Uo the phase should be 
close to zero. The median of the object phase is calculated 
only ones as y* n = median(angle(uQ)) without partitioning 
in two subsets as for the object amplitude. Finally, the filter- 
ing procedure by the BM3D filter (Steps 2 and 3 in SPAR - 
BC) is replaced by the following ones 

a^ +1/3 = BMZD a {abs{ul) - fa) + fa (35) 
a* +1/2 = BMZD a {*l +1/3 -p\) + p\ 
<^ +1/2 = BMZD v {angle{ul) - <p*J 

t+l t+l/2 / ■ / i+1/2 . t \\ 

The presented results of the object reconstruction from ex- 
perimental data are obtained using this modification BM3D 
filtering. 

5.4. Initialization for SPAR — BC: reconstruction with- 
out background compensation 

Here we consider the object reconstruction from the experi- 
mental data obtained by the approach originated in [15, 39]. 
The complex amplitude of the disturbed object is obtained by 
two various algorithm: by the mentioned AL algorithm [34] 
and by the successive phase-retrieval algorithm described in 
[40]. The second iterative algorithm is close to the circular 
wave reconstruction originated in [12, 13]: the calculated am- 
plitude is replaced by the square root of the given noisy in- 
tensity, keeping the calculated phase (the initial guess for the 
phase is zero). For simplicity, we refer to the latter algorithm 
as the Falldorf-Agour (FA) algorithm. 

In Fig. 3 we present the comparison of the reconstruction 
imaging of the disturbed object obtained from experimental 
data by these different methods. In the right column the es- 
timate of the disturbed object phase and amplitude computed 
by the FA algorithm [40] are illustrated (see Fig. 3(b) and 
Fig. 3(d), respectively). In the left column the reconstructed 
disturbed object phase and amplitude found by AL [34] are 
shown (see Fig. 3(a) and Fig. 3(c), respectively). These re- 
sults are shown for 25 iterations of the phase-retrieval algo- 
rithms. The artifacts which definitely should be addressed to 
the background are clearly seen in these images. It can be 
seen that the amplitude estimate by AL is significantly over- 
smoothed comparing with the result by FA. It is manifested in 
partial suppression of the diffraction artifacts on the geomet- 
rical elements with some degradation of a smooth surface as 
well. The phase reconstructions are not flat and have certain 
errors in the regions of the digits and geometrical figures in 




Fig. 3. Reconstructions of the "disturbed" object computed 
from experimental data (left column) by AL and (right col- 
umn) by FA [40]. In the top row the amplitude reconstruc- 
tions are presented: (a) by AL, (c) by FA. In the bottom row 
we demonstrate the phase estimates: (c) by AL, (d) by FA. 
The object reconstruction by AL (a[j from (a) and ip from 
(c)) is used for the initialization uf] = a[j o exp(j • ip®) in the 
SPAR - BC algorithm. 

the amplitude. The phase by AL has stronger degradation of 
the phase comparing with the phase calculated by [40], be- 
cause of a weak correction of the object reconstruction by 
Lagrange multipliers. We take quite small a = 1/20 because 
we are looking for a sharp result. It is found that larger a 
may correct the phase estimate but leads to more noisy am- 
plitude reconstruction, and denoising by larger regularization 
parameter fi (see Eq. (19)) results in oversmoothing. 

Note that the imperfect AL object estimate is used only as 
an initialization uj] = a§ o exp(j • <J? ) for the main procedure 
of SPAR - BC. 

5.5. SPAR — BC: reconstruction with background com- 
pensation 

The background reconstruction is produced with the calibra- 
tion experiments for the free space object u [k] = 1. Figure 
4(a) demonstrates the reconstructed background amplitude. 
However, the obtained results are appeared quite noisy and 
an additional postfiltering of the background amplitude a^ is 
introduced. Here we use the smoothed version of this back- 
ground reconstruction with no high frequency components. 
For the filtering of the amplitude and the phase of ub we 
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Fig. 4. Background compensation in SPAR — BC: (a) the 
reconstructed amplitude of the background s6s(ub) and (b) 
the initial guess for the object amplitude a6s(Ug) found with 
the smoothed background amplitude BA13D a (abs(uB))- 



again use BM3D filter [32]. The cross-sections of the original 
slb = afcs(us) and smoothed ag = BM3D a (a.B) back- 
ground amplitudes are illustrated in Fig. 5. The result of 
the compensation of the initial object amplitude by such an 
smoothed background a§[fc] = a&s(uQ[fc])/ae[fc] is shown 
in Fig. 4(b). The corresponding cross-section is presented in 
Fig. 8. 

The results of the wave field reconstruction obtained by 
the developed SPAR — BC algorithm are shown in Fig. 6 
for 25 iterations. In Figs. 6(a) and 6(b) the object ampli- 
tude reconstruction calculated with the smoothed and origi- 
nal background estimates are presented, respectively. It can 
be seen that the imaging of the amplitude estimate obtained 
using the smoothed background (Fig. 6(a)) is essentially bet- 
ter comparing with the reconstruction found with the back- 
ground without postfiltering, namely: the artifacts are well 
seen of the border of Fig. 6(b). Further improvement of the 
imaging can be achieved by postfiltering of the reconstructed 
Uq, and in Fig. 6(c) the result of such an additional BM3D 
filtering of u (computed again with the smoothed a^) is 
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Fig. 5. Cross-sections of (thin curve) the amplitude esti- 
mate of the background, as = a6s(us), and (thick curve) its 
smoothed version computed by BM3D filter, BM3D a (&B)- 
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Fig. 6. Object reconstruction by SPAR — BC. Compari- 
son of the reconstructed object amplitudes abs(u n ) calculated 
(a) with the smoothed background estimate and (b) with the 
original background reconstruction. The result of postfilter- 
ing of the object amplitude BM3D a (abs(uo)), r a j a = 0.04, 
is demonstrated in (c). The object phase estimate angle(uo) 
is illustrated in (d). The noise in the filtered phase estimate is 
totally suppressed. 

demonstrated. The reconstructed object phase angle(u.o) is 
illustrated in Fig. 6(d). 

The threshold parameter of the BM3D filtering in SPAR- 
BC is T a la = 0-13 for the object amplitude and r ¥ ,7„ = 2 
for the phase. Note that even with a large thresholding 
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Fig. 7. Cross-sections of the reconstructed object phase 
cp = angle(\io) and the filtered phase BM3D v (ip ), 
T <plip = 0-08. These results are shown along the dashed lines 
shown in Fig. 6(d). 
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Fig. 8. Cross-sections of (solid thin curve) the initial guess 
a.Q = <z&s(uq), (solid thick) the reconstructed amplitude 
So = abs(\io) and (dashed curve) the filtered amplitude esti- 
mate BM3D a (ko). These results are related to the imaging 
presented in Figs. 4(b), 6(a) and 6(c), respectively. The cross- 
sections are given along the dashed line in Fig. 6. 

T ¥>7y — ^ we have significant noise in the phase estimate (see 
Fig. 6(d)). Taking into account that ip [k] = 0, the recon- 
struction accuracy for the phase is RMSE=0.2. However, we 
can completely wipe the phase noise out by the mentioned ad- 
ditional filtering by BM3D with quite a small T^7„ — 0.08: 
compare the cross-sections of cp Q and BM3D v (ip ) in Fig. 
7. 

It is shown in Eq. (29) that the object reconstruction is a 
weighted sum of a noisy estimate from the propagation model 
and filtered approximation from the previous iteration. Thus, 
at each iteration the filtered object estimate is corrupted due to 
the used noisy intensity observations. It is a challenge to find 
a proper balance between denoising and oversmoothing. We 
obtain quite a sharp results of ao with some small distortions 
(see Fig. 6(a) and the related cross-section in Fig. 8). However 
remaining noise and artifacts can be additionally suppressed 
by BM3D filter after the main procedure of SPAR - BC. 
It provides crisp imaging, but the resulting amplitude is over- 
smoothed and small details are almost lost: the result of post- 
filtering of ao by 12 iterations with r a j a = 0.04 is presented 
in Fig. 6(c) and the corresponding cross-section in Fig. 8. 

6. DISCUSSION AND CONCLUSION 

It can be seen that the modified BM3D filtering (Eqs. (35)) 
works here as a classifier for the noisy binary object estimate. 
The estimate of the amplitude levels are found using the Otsu 
method, but the BM3D filtering shifts the value of the pixel 
ao[/c] to one of these two levels /3 or f3 x depending on the 
local neighborhood. 

The result of classification is presented in Fig. 9. Let 
if(-) stands for a histogram of a discrete distribution. The 
histogram if (ag) for the initial estimate of the compensated 
object amplitude is denoted in Fig. 9 by a solid curve. It is 
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Fig. 9. Partition produced according to the modified BM3D 
filtering (Eqs. (35)). (solid curve) if(ajj) is the histogram 
for the initial estimate of the object amplitude, a§ = a&s(ug) 
presented in Fig. 4(b). (dashed curve) if (ao) is the histogram 
for the resulting object amplitude estimate after 25 iterations, 
ao = abs(\io) illustrated in Fig. 6(a). Otsu's threshold with 
the upper and lower levels of the initial (p°, /3°, /3 ) and re- 
sulting object estimates (p, $ 1 , $ ) are also presented. 

already seen that the object ao looks to be binary. The his- 
togram of the resulting a is denoted here by a dashed curve. 
In Fig. 9 we also present the initial guess of the lower and 
upper levels f3\, the resulting lower f3 and upper levels 
and the Otsu threshold for the initialization p° and the final 
step p of our reconstruction. Note that blurred regions (see 
small details and borders on the geometrical elements e.g. in 
Fig. 6(a)) correspond to a "sloping valley" between two peaks 
of if (ao) in Fig. 9. 

In this paper a novel phase-retrieval algorithm with back- 
ground compensation and powerful BM3D filtering is pre- 
sented. The SPAR — BC algorithm demonstrates a very 
good reconstruction quality : we have a clear separation of the 
binary true object, and the background estimate "undertakes" 
strong fluctuations, which would be difficult to compensate 
by filtering only. The reconstructions by two different phase- 
retrieval methods (AL and FA) are presented to emphasize 
the obtained enhancement of imaging of the developed algo- 
rithm with respect to modern phase-retrieval algorithms with 
no background compensation (compare the results in Figs. 3 
and Fig. 6). 
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